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SECTION  1 
INTRODUCTION 


Sharp  pressure  waves  can  originate  in  aircraft  gas-turbine  engines  from  a 
variety  of  sources,  such  as  the  changing  inlet  conditions  due  to  a  sudden  maneuver, 
the  ingestion  of  ordnance  exhaust,  or  the  detonation  of  a  nearby  explosive.  These 
traveling  waves  interact  with  the  various  components  inside  the  engine  (compressor  and 
turbine  blade  rows,  burner  assemblies,  nozzles  and  ducts)  and  the  result  of  the  interaction 
is  to  produce  transmitted  and  reflected  waves.  Under  certain  circumstances,  the 
combined  effects  of  the  multiple  waves  produced  by  such  interaction  can  cause  major 
transients  in  the  engine  operation,  such  as  rotating  stall,  surge,  or  extreme  excursions 
in  the  combustor  temperature. 

Recent  measurements  (1-3)  provide  a  detailed,  time-resolved  data  base  for 
several  engines  under  carefully  controlled  conditions.  These  measurements  show  that 
a  pressure  wave  entering  the  engine  inlet  makes  its  way  through  the  various  stages  of 
the  machine  at  approximately  the  speed  of  sound  relative  to  the  local  flow  speed, 
provided  the  static  pressure  rise  is  on  the  order  of  ten  percent  or  less  of  the  initial 
static  pressure.  However,  the  amplifications  and  attenuations  of  the  incident  pressure 
pulse  that  occur  upon  passage  through  the  various  components  have  not  been  explained 
in  a  satisfactory  manner.  The  only  theoretical  tools  available  for  interpreting  the 
measurements  are  computer  codes  which  attempt  to  calculate  the  one-dimensional  time- 
dependent  flow  through  the  various  stages  (4-6).  None  of  the  available  codes  have 
been  applied  successfully  to  the  data  of  References  (1-3),  and  in  fact  most  of  them 
have  stability  and  waveform  restrictions  which  have  thus  far  prevented  their  application 
to  obtain  predictions  that  would  be  helpful  in  understanding  the  measurements. 

The  work  reported  here  has  made  some  progress  toward  the  prediction 
and  understanding  of  the  changes  in  waveform  and  amplitude  that  occur,  by  resorting 
to  an  analytic  solution  for  the  fundamental  problem,  namely,  each  compressor  or  turbine 
stage  is  represented  as  an  actuator  disk,  and  the  solution  is  found  for  the  waves  that 
are  transmitted  and  reflected  when  a  given  pressure  wave  is  incident  on  the  disk. 

The  report  is  organized  as  follows:  Section  2  contains  the  actuator-disk 
analysis  and  numerical  results  for  a  number  of  typical  cases.  Section  3  presents  the 
analytical  formulas  that  apply  for  very  weak  waves,  while  Section  k  treats  the  problem 
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of  multiple  reflections  between  a  pair  of  disks  as  a  means  of  estimating  how  an  actual 
engine  would  respond  to  a  disturbance.  Section  4  also  contains  a  comparison  with  some 
experimental  results,  which  are  more  easy  to  interpret  in  the  light  of  the  present  analysis. 

The  real-engine  experiments  include  a  number  of  features  that  are  not 
accounted  in  this  analysis,  and  the  Section  5  presents  a  discussion  of  these  features 
and  the  'jrther  research  that  is  now  made  possible  by  the  results  of  this  work. 
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SECTION  2 

ACTUATOR-DISK  ANALYSIS 


The  data  presented  in  (1-3)  include  detailed  pressure  measurements  at  a 
variety  of  points  inside  operating  engines  through  which  shock  waves  have  been  sent. 
Pressure  histories  were  recorded  on  gages  placed  in  positions  upstream  and  downstream 
of  various  engine  components:  the  inlet  fan,  the  core  compressor,  etc.  From  these 
data  it  is  possible  to  infer  how  an  incident  shock  wave  interacts  with  each  component; 
the  pressure- time  records  at  the  various  positions  reveal  a  very  complex  system  of 
compression  and  expansion  waves.  Attempts  to  interpret  these  complex  results  have 
been  thwarted  in  the  past  because  of  the  absence  of  a  solution  for  the  elemental 
problem  that  arises  in  such  a  flow,  namely  the  interaction  of  the  shock  with  a  single 
component.  This  fundamental  problem  is  the  one  whose  solution  is  reported  herein. 

2.1  ANALYSIS. 

In  order  to  make  the  problem  manageable,  each  comjjonent  is  treated  as 
an  actuator  disk  (7),  i.e.  the  axial  extent  of  the  component  is  assumed  to  be  shortened 
into  a  sheet  of  zero  thickness,  across  whidi  the  initial  and  final  states  are  connectd  by 
a  discontinuity.  This  approximation  is  suitable  as  long  as  one  examines  the  solution  on 
a  time  scale  that  is  large  compared  to  the  time  for  a  sound  wave  to  pass  through  the 
stage.  During  the  latter  time,  details  of  the  inside  of  the  stage  are  important;  for 
times  large  compared  to  this,  however,  the  action  of  the  stage  is  adequately  represented 
by  its  inlet  and  outlet  conditions,  considered  to  take  place  by  a  discontinuous  jump  in 
a  region  of  vanishing  axial  extent. 

The  representation  of  a  turbomachine  stage  by  a  disk  of  zero  axial  extent 
across  which  discontinous  changes  in  stagnation  conditions  occur  is  an  excellent 
approximation  whenever  details  of  the  flow  inside  the  stage  are  not  important.  This 
approximation  has  been  applied  with  great  success  to  a  number  of  problems.  For 
example,  the  propagation  of  plane  waves  through  ducts  has  been  studied  for  many  years, 
and  solutions  are  available  for  the  changes  that  occur  when  such  waves  encounter 
discontinuities  of  various  kinds,  such  as  sudden  area  changes,  flame  fronts,  and  other 
moving  waves  (8).  Surprisingly,  the  solution  for  the  case  where  the  discontinuity  is  an 
actuator  disk  appears  not  to  have  been  studied,  although  the  basic  equations  applicable 
to  this  case  are  discussed  in  several  papers  (9-11). 
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2.1.1  Qualitative  Results. 

The  present  paper  contains  the  solution  for  the  waves  that  are  generated 
when  a  pressure  wave  (either  a  shock  or  an  expansion)  arrives  at  an  actuator  disk. 
The  disk  is  characterized  as  having  a  given  stagnation  pressure  ratio,  and  is  assumed 
to  produce  this  change  isentropically.  The  types  of  the  reflected  and  transmitted  waves 
are  found  (i.e.,  whether  they  are  expansions  or  compressions),  in  addition  to  the 
magnitudes  of  the  changes  they  produce  relative  to  that  of  the  incident  wave.  It  is 
found  that  a  shock  wave  incident  from  upstream  on  a  compressor  stage  reflects  as  an 
expansion  and  is  transmitted  as  a  weakened  shock.  An  expansion  wave  incident  from 
upstream  reflects  from  a  compressor  stage  as  a  shock,  and  is  transmitted  as  a  weakened 
expansion. 


These  results  are  reversed  when  the  wave  is  incident  from  the  downstream 
side  of  a  compressor  stage:  an  incident  shock  reflects  and  is  transmitted  as  a  shock 
while  an  incident  expansion  produces  both  a  transmitted  and  reflected  expansion.  These 
results  are  in  agreement,  both  qualitatively  and  quantitatively,  with  the  measurements 
of  (1-3).  The  sketch  below  contains  a  qualitative  summary  of  the  transmission  and 
reflection  laws.  (There  are  eight  cases,  corresponding  to  two  values  of  the  stagnation 
pressure  ratio  (compressor  or  turbine),  two  signs  of  the  incident  pressure  rise  (shock 
or  expansion  wave),  and  two  directions  in  which  the  incident  wave  moves  (downstream- 
moving  or  upstream-moving)).  The  sketch  contains  the  wave  diagrams  for  the  various 
cases.  The  diagrams  are  shown  in  a  graph  whose  vertical  axis  is  the  time  and  whose 
horizontal  axis  is  axial  position,  in  accordance  with  the  conventional  methods  of  one¬ 
dimensional  gas  dynamics  (8). 


COMPRESSORS 


INCIDENT  SHOCK 


TURBINES 


INCIDENT  SHOCK  '  INCIDENT  EXPANSION 


surge  and  rotating  stall  (10,  11).  Those  developments  pertain  generally  to  engine 
transients  of  a  much  longer  time  scale  (the  order  of  the  engine  length  divided  by  the 
speed  of  sound),  and  make  use  of  unsteady  compressor  characteristics  in  the  form  of 
instantaneous  pressure  rise  vs.  mass  flow  relations.  The  present  analysis  is  capable  of 
being  extended  to  incorporate  similar  characteristics,  but  now  for  individual  stages 
rather  than  for  the  whole  component.  Such  an  application  would  assist  in  showing  how 
surge  and  rotating  stall  develop  within  the  comjsonent. 

These  analytic  results  can  also  be  used  to  simplify  the  computer  codes 
referred  to  above,  and  should  be  capable  of  extending  their  range  of  applicability.  This 
is  because  the  codes  have  difficulty  in  tracking  the  large  numbers  of  waves  that  are 
generated;  the  ability  to  describe  the  locations  and  amplitudes  of  these  waves  analytically 
may  enable  significant  improvements  in  the  numerical  modeling  of  these  events. 

2.1.2  Solution  Method 

Once  the  actuator-disk  approximation  is  made,  the  solution  of  the  problem 
can  be  carried  out  directly  using  the  classic  techniques  of  wave-diagram  construction 
(8).  The  Riemann  invariants  which  form  the  basis  of  the  solution  are  the  same  as 
those  for  flow  in  a  duct,  and  need  only  be  modified  for  the  transitions  across  the 
actuator  sheet.  If  one  were  interested  in  details  of  the  flow  within  an  engine  component, 
the  analysis  would  have  to  be  generalized  to  include  the  effects  of  angular  velocity 
changes  within  the  component,  which  would  then  be  referred  to  as  an  'actuator  duct'. 
A  brief  discussion  of  these  levels  of  approximation  can  be  found  in  (7)  and  (9). 

For  the  actuator-disk,  one-dimensional-duct-flow  model,  the  flow  can  be 
charted  in  an  axial  distance  (x)  vs  time  (t)  diagram  (see  Figure  1).  (Here  the  flow 
quantities  are  constant  in  various  zones,  which  are  separated  by  moving  waves.)  Thus, 
the  vertical  axis  represents  the  actuator  disk  and  the  time  axis  and  the  lines  separating 
regions  1  and  2  or  4  and  5  are  the  incident  and  transmitted  shock  waves.  The  broken 
line  between  regions  5  and  6  is  the  interface  which  separates  gas  that  has  or  has  not 
been  affected  by  the  reflected  wave  that  separates  regions  2  and  3.  The  reflected  wave 
is  shown  as  a  shock  in  Figure  1;  the  results  below  will  show  that  it  is  actually  an 
expansion  wave.  Particle  paths  are  shown  as  the  short  dashed  lines. 
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Figure  1.  Wave  diagram  for  a  dowrtstream-propagating  wave. 


The  actuator  disk  is  assumed  to  have  a  given  constant  stagnation-pressure 
ratio  across  it,  and  to  achieve  its  compression  isentropically,  i.e., 


T 

r  =  ^  =  (n)  Y 


(1) 


where  subscripts  A  and  ©  denote  conditions  ahead  of  and  behind  the  disk,  Pqj  and  Toi 
are  the  stagnation  pressure  and  temperature  in  region  i,  iT  is  the  ratio  of  specific 
heats  (assumed  constant),  and  ff  is  a  given  constant.  A  second  condition  across  the 
actuator  disk  is  that  the  mass  flow  be  constant,  i.e. 
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(2) 

where  P  and  U  denote  the  density  and  axial  velocity,  respectively. 

The  solution  procedure  is  as  follows: 

1)  The  compression  ratio  Tf  and  specific-heat  ratio  if  are  chosen;  the 
initial  inlet  Mach  number  Mj  is  chosen  and  the  strength  of  the  incident  wave  (as  given 
by  the  pressure  ratio  P2/P1)  chosen.  Conditions  across  the  incident  wave  are  given 
by  either  the  shock-wave  or  expansion-wave  relations.  For  the  shock-wave  case,  these 
are: 
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where  K  5^^  denotes  the  velocity  of  the  shock  wave  connecting  regions  A  and  B.  The 
first  of  these  gives  the  shock  velocity,  the  second  the  particle  velocity. 


7 


For  the  case  of  an  expansion  wave  incident  from  the  upstream  side,  the 
conditions  behind  the  wave  are  found  from  the  constancy  of  the  Riemann  invariant  Q 
across  the  wave; 


Q  -  ^  ^5  ~  ^3  - 


(7) 


where  the  sound  speed  Ou  can  be  found  from  the  isentropic  relation 

(Lb 


zt 


a 


The  remaining  equations  are  then; 


(8) 
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The  leading  edge  of  the  expansion  wave  travels  at  the  speed: 

dx/dz  =  1/^  +  ^2^ 


while  its  trailing  edge  travels  at  the  speed  t/lg  +  CLq 


2)  Conditions  behind  the  actuator  disk  are  found  from 


-1 


(10) 


(11) 


^'P  '^5  T  ~ 


r  T  ^  ^ 


(12) 


By  using  the  isentropic  relation: 
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T'bAa  =  (ps/PaA  '  = 


(13) 


these  can  be  put  into  the  single  equation 


This  equation  can  be  solved  by  successive  substitutions,  using  =  0  as  an  initial 

guess. 


It  should  be  noted  that  the  designations  "ahead  of"  and  "behind"  are  to 
be  interpreted  for  the  actuator  disk  as  referring  to  a  pair  of  points  which  lie  on  the 
upstream  and  downstream  sides  of  the  disk,  at  the  same  instant  of  time  on  the  wave 
diagram.  When  these  designations  are  applied  to  moving  shocks  and  expansions,  they 
are  to  be  interpreted  as  denoting  points  which  are  at  a  fixed  spatial  position  and  are, 
respectively,  at  times  less  than  and  greater  than  the  time  at  which  the  wave  passes 
the  given  position. 

3)  Next,  values  of  P5/P4  and  P3/P2  are  chosen,  and  these  are  adjusted 
(by  an  iteration  process  described  below)  until  values  are  found  such  that  pg  =  p^ 
and  Kg  =  U5  •  The  relations  between  conditions  in  regions  5  and  4  are  given  by  either 

the  shock-  or  expansion-wave  relations  above,  with  A  =  4  and  B=5. 

The  transition  2-3  may  be  a  reflected  shock,  in  which  case  the  shock 

relations  are  used,  with  A=2,  B=3.  Because  this  wave  is  moving  to  the  left,  care  must 

be  taken  to  choose  the  appropriate  signs;  the  positive  square  roots  are  taken,  and  the 
equations  for  the  shock  and  particle  velocity  are  written  as: 
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In  both  of  these  expressions  the  quantity  t/g  is  a  negative  number. 

If  the  transition  2-3  is  a  (left-running)  expansion  wave,  the  solution  is 
given  by  the  constancy  of  the  Riemann  invariant  P  : 

III 

This  can  be  solved  explicitly  as 


The  speed  of  this  expansion  wave  varies  between  the  limits; 

=  ^6-^6  (19) 

d-t  ^  ^  dt 

The  iteration  process  used  to  update  the  values  of  p^/pit  and  P3/P2  ‘s 
based  on  a  straight-line  extrapolation;  the  graphs  of  p^  vs  u^  and  pg  vs  ug  found  for 
various  values  of  p^/p^  and  P3/P2  are  nearly  straight  lines.  Thus,  their  point  of 
intersection  can  be  readily  calculated,  and  the  corresponding  values  of  pj/p^  and  P3/P2 
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used  for  the  subsequent  iteration.  These  iterations  were  terminated  when  successive 
approximations  differed  by  less  than  10“^. 

For  cases  where  the  incident  wave  comes  from  downstream,  a  few 
modifications  to  these  relations  are  required.  Figure  2  shows  the  wave  diagram,  which 
now  has  seven  zones.  As  before,  M|,  Tf  and  ^  are  given,  and  the  relation  across 
the  actuator  disk  between  regions  1  and  4  is  given  by  equation  14,  with  A=l,  B=4. 


K 


Figure  2.  Wave  diagram  for  an  upstream-propagating  wave. 

The  incident  pressure  ratio  p^/p^j  is  also  given;  if  it  is  greater  than  1.0, 
equations  3  and  4  are  used,  corresponding  to  an  upstream-moving  shock;  if  it  is  less 
than  1.0,  equations  15  and  16  are  used,  for  an  upstream-moving  expansion  wave. 

The  transition  5-6  uses  the  relations  for  a  right-moving  shock  (equations  3- 
4)  or  expansion  (equations  7-8),  while  that  between  regions  1  and  2  uses  the  relations  for 
a  left-moving  shock  (equations  15-16)  or  expansion  (equation  17).  Regions  2  and  7  are 
connected  by  the  actuator-disk  relation,  equation  4,  with  A=2,  B=7. 

An  iteration  process  similar  to  that  used  for  waves  incident  from  the 
upstream  side  is  now  used:  a  series  of  wave  strengths  pg/p^  produces  a  graph  of  ug 
versus  pg,  while  a  series  of  strengths  P2/P1  produces  a  graph  of  Uy  versus  py.  Linear 
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fits  to  these  curves  are  used  to  find  improved  values  of  pg  and  P2  (the  latter  requiring 
a  further  linear  fit  to  the  variation  of  P2  versus  p^). 

Fortran  programs  incorporating  these  formulas  are  given  in  Appendices  A 

and  B. 


2.2  NUMERICAL  RESULTS. 

Consider  first  the  case  of  a  shock  wave  incident  from  the  upstream  side  of 
a  compressor  stage.  Tables  1-4  present  typical  results:  the  first  one  shown  is  a 
baseline  case,  chosen  for  its  similarity  to  the  first- stage- fan  data  of  (1).  The  subsequent 
tables  then  show  the  effects  of  changing,  one  at  a  time,  the  parameters,  Tf,  M  j  and  P2/pi* 


Table  1.  Numerical  solution  for  77  =  2,  *  0.2,  ^  =  1.4,  P2/Pi  *  1.1 


i  s 

2 

J 

6 

5 

6 

p./pi 

l.OW* 

2.0397 

2.1950 

2.1951 

J  i 

1.070'. 

1.0350 

1.6639 

1.7536 

1.7536 

T,/Ti 

1.0276 

1.0139 

1.2259 

1.2519 

1.2519 

U,/3l 

.2636 

.3025 

.1202 

.1786 

.1786 

M, 

.2669 

.3006 

.1086 

.1596 

.1596 

a, /a, 

1.0137 

1.0069 

1.1072 

1.1189 

1.1189 

.9999 

.9999 

1.9999 

1.9998 

"  1.566 
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Table  2. 

Numerical  solution  for 

“4.M, 

=  0.2, 

1.4, 

P2/P1  *  1.1 

‘  ' 

2 

3 

4 

5 

6 

Pi/pi 

1.1 

1.0061 

4.1029 

4.3170 

4.3170 

J  1 

1.0704 

1.0043 

2.7411 

2.8425 

2.8424 

T,/Tl 

1.0276 

1.001$ 

1.4968 

1.5187 

1.5188 

u,/a  1 

.26S6 

.3327 

.0730 

.1176 

.1176 

.2649 

.3324 

.0596 

.0954 

.0934 

a, /a, 

1 .0 1 37 

1.0009 

1.2234 

1.2324 

1.2324 

Poi^Po  1 

.9999 

.9999 

4.0 

3.9999 

3.9996 

^^'3.6 

(eu),,2 

■■  1.671 

Table  3. 

Numerical  solution  for  TI 

=  2.Mi 

-0.2,V= 

1.4, 

P2/P1  =  1.5 

X  5 

2 

3 

4 

5 

6 

p/pi 

1.3 

.9914 

2.0397 

2.8149 

2.8149 

J  1 

1.3333 

.9919 

1.6639 

2.0923 

2.0923 

T,/T, 

1.1250 

.9995 

1.2259 

1.3454 

1.3467 

u,/j  J 

.49SS 

.$035 

.1202 

.3813 

.3813 

.4703 

.8037 

.1086 

.3287 

.3285 

a,, 'a; 

1.0607 

.9997 

1.1072 

1.1399 

1.1603 

Poi'^Po) 

.9932 

.9932 

2.0 

1.9932 

1.9865 

(9u),,2 

'  3.933 

Table  4. 

Numerical  solution  for  TI 

=  2.M, 

=  0.4,  = 

1.4 

,  P2/Pi  =  1.1 

1  : 

2 

3 

4 

5 

6 

Pi/Pi 

1.1 

1.012$ 

2.1680 

2.3374 

2.3374 

i/  1 

1.0704 

1.0091 

1.7380 

1.8339 

1.8339 

T,/T, 

1.0276 

1.0037 

1.2474 

1.2746 

1.27  46 

Ui/3, 

.4686 

.5280 

.2302 

.29  0  5 

.29  0  5 

M| 

.4622 

.5270 

.2061 

.2573 

.2573 

ai/a, 

1.0137 

1.0018 

1.1169 

1.1290 

1.1290 

Poi/Pol 

.9999 

.9999 

2.0 

1.9999 

1.9998 

"  1.332 

Cu)|_2 
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Several  items  are  worthy  of  note: 


1)  The  reflected  wave  is  an  expansion.  It  might  have  been  expected  that  the  wave 
reflected  from  a  high-solidity  compressor  face  would  be  a  shock  wave,  but  these  results 
show  clearly  that  it  is  an  expansion  wave,  for  the  whole  range  of  cases  studied.  This 
finding  is  also  confirmed  by  the  pressure/time  traces  of  (I). 

2)  The  "transmission  coefficient",  defined  as 

T  -  (20) 

^  -  -Pz-p, 

lies  in  the  range  0.5  to  0.8  of  ff  ,  favoring  the  0.8  Tf  value  when  ff  itself  is  2.0. 
This  result  compares  favorably  with  the  values  presented  in  (1)  (see  Figs.  19  and  20). 

3)  The  arrival  of  the  incident  shock  causes  a  very  large  change  in  mass  flow  rate  -  on 
the  order  of  a  factor  of  1.5  -  despite  the  fact  that  the  incident  wave  is  a  relatively 
weak  shock. 

4)  Other  calculations,  not  shown  here,  reveal  that  these  qualitative  conclusions  are 
not  sensitive  to  variations  in  the  specific-heat  ratio. 

The  agreement  between  the  predicted  and  observed  values  of  the 
transmission  coefficient  is  encouraging,  while  the  indication  that  the  reflected  wave  is 
an  expansion  is  in  qualitative  agreement  with  the  results  of  (1)  and  (3).  The  conditions 
of  the  latter  experiment  contain  other  factors  that  affect  the  pressure  signals;  these 
are  discussed  at  greater  length  below.  Overall,  these  results  are  a  useful  guide  in 
interpreting  the  response  of  individual  stages. 

The  large  mass  flow  perturbation  raises  the  question  whether  the  actuator- 
disk  stagnation-pressure  ratio  can  properly  be  regarded  as  constant.  While  it  is  true 
that  the  engine  RPM  does  not  change  on  the  time  scale  of  these  wave  transits, 
nevertheless  the  unsteady  compressor  characteristics  used,  for  example,  in  the  work  of 
Moore  and  Greit2er(12, 13)  indicate  that  the  quantity  rf  should  be  allowed  to  vary 
with  P  3U3.  Very  little  is  known  at  present  about  the  variation  of  pressure  with 
mass  flow  on  the  relatively  short  time  scale  of  a  shock  passage  through  a  stage,  and  it 

is  not  clear  that  the  longer-time  variations  used  in  the  rotating-stall  work  cited  above 
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are  accurate  for  the  present  problem.  Further  investigation  of  the  detailed  mass- 
flow/pressure  variations  during  wave  passage  through  a  single  stage  is  needed  to  resolve 
this  question. 


Results  of  further  calculations  are  given  in  Figures  3-4  for  shocks  and 
expansions  which  are  incident  both  from  upstream  and  downstream.  These  figures 
confirm  the  general  transmission  and  reflection  laws  mentioned  in  the  Introduction. 


®/ 

© 

Figure  3.  Actuator-disk  solutions  forx  =  2.0,  M-  =  0.2,  7  *  1.4 
(Downstream -propagating  waves). 


Figure  4.  Actuator-disk  solutions  for  x  *  2.0,  M-j  =  0.2,  7  *  1.4 
(Upstream-propagating  waves). 


Similar  results  have  been  found  for  the  transmission  and  reflection  of 
sound  waves  in  turbomachinery  components  (14-15).  The  latter  studies  are  somewhat 
more  complex,  in  that  they  consider  sinusoidal  waves  and  allow  for  wavefront  propagation 
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in  directions  other  than  parallel  to  the  axis  of  rotation.  For  the  case  of  axial 
projiagation,  these  papers  reach  the  same  conclusions  about  transmission  and  reflection 
as  those  derived  above.  The  present  results  show  that  there  is  no  qualitative  change 
in  the  transmission  and  reflection  laws  for  finite  wave  amplitudes. 

There  are  limits  to  the  validity  of  the  present  solution,  however.  The 
curves  in  Figures  3-4  stop  at  finite  limits;  for  example,  in  the  case  of  waves  incident 
from  the  upstream  side,  a  sufficiently  strong  reflected  shock  will  stop  the  flow  (M3  =  0), 
while  a  sufficiently  strong  reflected  expansion  wave  will  choke  the  flow  (M3  =  1). 
Beyond  these  limits,  there  is  no  solution  of  the  type  being  considered  here.  It  is 
interesting  to  note  that  these  limits  are  reached  for  fairly  moderate  values  of  the 
incident  wave  strength;  for  example,  an  incident  expansion  wave  with  a  pressure  drop 
of  as  little  as  twenty  percent  of  the  initi2d  pressure  is  sufficient  to  stop  the  flow. 
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SECTION  3 

ANALYTICAL  SOLUTIONS  FOR  WEAK  WAVES 


The  incident-wave  strengths  of  interest  in  many  applications  have  pressure 
ratios  close  to  one;  thus  it  is  useful  to  examine  the  analytic  forms  that  the  full 
equations  reduce  to  in  this  limit. 


Consider  first  the  case  where  the  initial  disturbance  is  incident  from 
upstream:  >  ^^1>  ^  given.  Let  the  (given)  incident  wave  strength 

be  denoted  by 


1  +  £ 


£  1 


(21) 


It  can  be  expected  that  the  solution  will  have  the  form 


;  6  (£)«  ^  ^  <sc  ^ 


(22) 


Where  the  dependence  of  6 
^6  =  ^5  >  ^6  ^  ^5  . 


and  y^on  €  is  to  be  found  from  the  usual  matching  condition 


The  quantities  £  ,  S>  ,  and  "V  can  all  be  positive  or  negative,  with  no 
change  in  the  functional  form  of  the  equations,  since  weak  shocks  and  expansions  are 
indistinguishable  to  the  third  order  in  small  disturbances. 


The  desired  equations  are  found  by  expanding  the  exact  relations  in  Taylor 
series  (with  the  exception  of  the  transition  across  the  actuator  disk  -  see  below).  In 
these  expansions,  the  Mach  numbers  are  not  assumed  to  be  small,  adthough  their 
numerical  contributions  will  in  some  cases  be  negligible. 

The  pressures  in  regions  5  and  6,  which  are  required  to  be  equal,  are 

(23) 
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=>5  P4 


7  =  77  ('"V 


-  J1  ^  IL  -  ^ 

■Pi  -^i  Pi  ^ 


t  +  (S'  ^  + 


--) 


(24) 


Thus 


The  particle  velocities  in  regions  5  and  6  can  be  written  as: 


(25) 


^1 


Us  U4 


U4  Ui  af 


=  777 


U4  U5 


U4. 


U. 


=  r^( 


u 


1 


'V 


Vh. 


(26) 


a, 


u- 


Ux  a. 


u 

"Ur 


6 


£-S 


(27) 


Equating  these  two  gives 


Ui 


1  -h 


'V 


u, 


u- 


£-6 


(28) 


A  simultaneous  solution  of  equations  25  and  28  gives: 
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^1  VP5  'pj  ^<^4  ^  [w,  P,  U3  'VH^  Pg  U,_ 

1  ^  ^  ^  _L.  ^ 

VMj  P/  U3  ^4  P5  U| 


^  ^6  /^fe  ^4  \  '_L  ^  ^  %  ^6 

V  - - -'(30) 

y  ^6  /  ^6  ^4 

VM,  P;  tl3  ■'■  ifM4  P3  Uf 

These  equations  cannot  be  solved  directly,  because  the  pressure  and  velocity  ratios 
across  regions  3  and  6  require  knowledge  of  M3,  which  in  turn  depends  on  ^  and  6  : 


In  addition,  the  fact  that  the  actuator-disk  transition  requires  the  solution  of  an  implicit 
equation  prevents  the  derivation  of  an  explicit  solution  valid  for  arbitrary  stagnation 
pressure  ratios. 

However,  for  the  case  where  Mj  is  small  (and  of  the  same  order  as  €  ) 
a  fully  analytic  solution  is  possible.  In  this  limit. 


^1-1^4  _ 

M  f  +  M4  ^  ^  +  1^4 


(32) 


Moreover,  Mj  and  M/^  are  related  in  this  limit,  so  that  the  solution  takes  the  simple  form 


M4  =  n 


iv 


_  ^  _  y-»-/ 

1-Tf  _  2  77 

-111  ^  ^  ~  _ H[ 

1-^n  HTT  ^,'4 


(33) 


These  relations  contain,  qualitatively,  the  results  of  the  previous  section.  For  a 
compressor,  for  example,  is  less  than  Mj;  thus  shocks  reflect  as  expansions  and 
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vice  versa,  while  the  transmitted  wave  is  of  the  same  type  as  the  incident  wave,  with 
reduced  amplitude.  For  a  turbine,  is  greater  than  thus  incident  shocks  and 
expansions  are  reflected  as  shocks  and  expansions  of  reduced  amplitude,  and  are 
transmitted  as  shocks  and  expansion  of  increased  amplitude. 


The  derivation  of  these  results  is  based  on  approximating  the  solution  of 
the  (implicit)  actuator-disk  equation  for  small  values  of  the  upstream  Mach  number: 


Thus 


U 


B 


U 


=  r 


=  r 


=  t 


f 

x-1 


M  Z 

- 

+  0 


u 


B 


U 


(34) 


(35) 


By  substituting  this  back  on  the  right-hand  side  of  the  original  equation,  it  is  found  that 


(36) 


For  a  compressor,  this  equation  is  valid  within  three  percent  over  the  range 
for  y  =  1.4  and  Tt  =  2  (see  Figure  5).  The  approximation  is  not  as  accurate  for  a 
turbine,  where  Ug  »  fhe  results  in  Figure  6  show  that  it  is  accurate  to  a  few 

percent  over  about  two  thirds  of  the  allowable  range  of  M^.  This  corresponds  to  the 
range  for  which  the  Mach  number  downstream  of  the  disk  stays  less  than  about  0.75. 

This  approximate  solution  can  now  be  used  to  simplify  the  small-disturbance 
solution  (eqs.  29  and  30).  The  pressure  ratio  is  given  by 


E  ^  1  -h  r  y-i=/-n  ^ 
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actuator-disk  solutions 


-  EXACT 

- APPROXIMATE 


Figure  5.  Actuator-disk  solution  for  x  =  2.0.  Figure  6.  Actuator-disk  solution  for  x  -  0.6. 


When  eqs.  29  and  30  are  expanded  using  these  equations,  it  is  found  that  all  of  the 
pressure-  and  velocity-ratio  terms  cancel,  leaving  a  remainder  of  order  M|f  ,  which  is 
of  higher  order  than  the  terms  retained.  A  number  of  intermediate  solutions  can  also 
be  found,  if  certain  combinations  of  terms  are  retained  and  the  relative  orders  of  M| 
and  S  are  ignored.  Given  the  simplicity  of  the  numerical  solutions  themselves,  it  is 
recommended  that  they  be  used  where  accuracy  is  considered  important,  and  that  the 
analytic  solutions  be  used  where  qualitative  information  is  the  primary  requirement. 

For  the  case  where  the  incident  wave  comes  from  the  downstream  side, 
a  similar  analysis  leads  to 
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(38) 


P  = 

A  = 


M4  - 

+  M4 

M|  +  M4. 


where  _ 

P4 


1  ■+■  oc 


1-hM^  oc  give  a 


These  approximate  formulas  (eqs.  32  and  38)  confirm  the  transmission  and 
reflection  laws  summarized  in  Section  2. 
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SECTION  4 

MULTIPLE  REFLECTIONS  BY  A  PAIR  OF  ACTUATOR  DISKS 

As  will  be  seen  below,  the  above  analysis  is  in  general  agreement  with 
the  experiments  of  (1-3),  in  that  the  pressure  measured  upstream  of  the  compressor 
shows  a  reflected  expansion,  follovwng  the  passage  of  an  incident  shock.  The  detailed 
waveform  of  the  reflection  is  not  a  step  function  for  several  reasons,  among  them  the 
fact  that  the  incident  wave  is  not  a  step  function,  as  well  as  the  fact  that  the  reflected 
waveform  is  affected  by  multiple  reflections  from  the  various  stages  of  the  compressor. 

In  order  to  shed  some  light  on  the  latter  process,  calculations  have  been 
done  for  two  pairs  of  actuator  disks:  the  first  pair  consists  of  two  stages  of  compression, 
while  the  second  consists  of  a  compressor  followed  by  a  turbine. 

For  the  first  pair,  both  disks  were  taken  to  have  stagnation-pressure  ratios 
of  2.0  and  an  incident  shock  strength  p2/pi  =  1»2  was  used  with  an  initial  flow  Mach 
number  of  0.2  and  specific-heat  ratio  if  =  1.4  The  main  results  are  shown  in  Fig.  7, 
which  contains  pressure  and  mass- flow  data  in  the  various  regions  of  the  x,t  diagram. 

Before  discussing  these  results,  two  observations  are  necessary;  the  first 
is  that  all  entropy  discontinuities  have  been  ignored.  In  general,  the  transmitted  waves 
are  followed  by  interfaces  across  wWch  the  pressure  and  velocity  are  matched,  but 
there  is  a  discontinuity  in  entropy  (or,  equivalently,  in  local  stagnation  conditions)  across 
the  interface.  Because  the  wave  amplitudes  and  flow  Mach  numbers  studied  in  these 
cases  are  small,  it  is  legitimate  to  neglect  the  waves  that  would  normally  be  generated 
by  interactions  with  the  interface,  and  the  validity  of  this  assumption  is  borne  out  by 
the  fact  that  the  stagnation  quantities  calculated  here  are  essentially  constant  throughout 
the  region  of  multiple  reflections. 

A  second  observation  is  that  the  wave  speeds  have  all  been  taken  as  equal 
to  +a  in  Figures  7  to  9  to  simplify  the  presentation.  The  actuad  wave  speeds  vary  by 
+  20%  about  this  baseline  value,  and  expansion  waves  diverge  as  they  propagate,  but 
these  details  are  suppressed  in  the  presentation. 
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Figure  7.  Pressure  and  massflow  values  during  wave  reflections  by  two  compressor 
actuator  disks;  =  0.2,  7  =1.4. 


The  pressures  shown  in  Figure  7  for  the  two  stage  compressor  case  indicate 
that  the  inlet  pressure  returns  to  its  original  level  after  about  four  reflected  waves, 
and  undergoes  a  slight  undershoot  during  its  return.  The  mass-flow  data  show  a 
monotonic  increase  to  the  post-disturbance  value  on  the  upstream  side  of  the  compressor, 
while  the  values  downstream  overshoot  at  first  and  then  subside.  After  about  four 
reflections,  the  mass  flow  rate  has  been  increased  by  a  factor  of  2.36.  This  is  the 
same  result  (to  three  significant  figures)  as  would  be  obtained  by  a  single  actuator  disk 
of  stagnation  pressure  ratio  equal  to  4.0.  The  fact  that  the  mass  flow  can  be  sustained 
at  this  level  even  while  the  inlet  pressure  has  returned  to  its  pre- disturbance  value 
assumes  that  the  conditions  far  upstream  remain  unchanged;  i.e.,  no  waves  come  from 
upstream  after  the  initial  disturbance. 

Figure  8  displays  the  pressure-time  history  that  would  be  seen  at  a  distance 
L  ahead  of  the  first  disk.  It  bears  a  qualitative  resemblance  to  the  measurements  of 
(3),  and  serves  to  emphasize  that  the  reflected  system  corresponding  to  a  step  input 
will  be  a  series  of  expansions  and  compressions.  For  a  multistage  compressor,  this 
series  will  appear  to  be  nearly  a  continuous  variation. 
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Figure  8.  Static-pressure  history  at  a  station  upstream  of  a  pair  of  compressive  actuator  disks. 

Results  for  the  compressor/turbine  pair  of  actuator  disks  are  shown  in 
Figure  9,  using  the  same  format  as  above.  In  this  case,  there  is  an  overshoot  in  the 
inlet  mass  flow,  while  the  inlet  pressure  falls  at  first  and  then  returns  to  the  value 
generated  by  the  first  incident  wave.  These  final  values  are  achieved  after  about  four 
reflections;  further  waves  have  pressure  ratios  of  one.  The  return  of  the  pressure  to 
the  value  generated  by  the  first  wave  is  consistent  with  what  would  be  predicted  by 
collapsing  the  two  actuator  disks  into  a  single  one:  since  the  product  of  their  stagnation 
pressure  ratios  is  1.0,  it  is  as  though  no  disks  at  aUl  were  present,  and  the  long-term 
solution  remains  at  the  values  set  by  the  incident  wave. 
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Figure  9.  Pressure  and  massflow  values  during  wave  reflections  by  a  contpressor  turbine 
actuator-disk  pair;  M-j  =  0.2,  7=1.4. 

COMPARISON  WITH  EXPERIMENT. 


As  noted  earlier,  there  are  many  features  of  the  full-scale  experiment 
(1-3)  that  are  not  accounted  for  in  the  present  analysis.  Nevertheless,  several  of  the 
qualitative  features  of  the  measurements  are  the  same  as  those  noted  here.  Figure 
10  shows  the  pressure- time  histories  recorded  at  several  stations  upstream  of  the  inlet 
of  the  TF33  engine.  These  data  are  replotted  in  Figure  II  so  as  to  show  the  pressure 
waveforms,  at  three  instants  of  time.  Note  first  that  at  10  milliseconds  (zero  time 
corresponds  to  the  instant  when  the  incident  wave  emerges  from  the  shock  tube)  the 
incident  wave  is  not  a  step  function,  but  rather  shows  an  increasing  amplitude  with 
distance  from  the  engine  inlet  (x  =  0).  At  t  =  20  milliseconds,  however,  the  pressure 
nearest  the  engine  inlet  has  dropped,  and  this  lower-pressure  region  appears  to  have 
moved  through  the  incident  pulse  on  the  30-millisecond  trace. 

The  fact  that  the  incident  wave  is  not  a  step  function  is  due  to  several 
effects;  one  of  the  primary  ones  is  the  area  change  at  the  shock-tube  exit  euid  also  that 
at  the  engine  inlet.  It  is  virtually  impossible  to  maintain  a  strictly  flat  pressure  profile 
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Figure  10.  Pressure  time  histories. 
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in  this  experiment,  and  a  rigorous  interpretation  of  the  data  requires  a  more  detailed 
calculation  in  which  the  evolution  of  the  incident  waveform  can  be  followed.  Moreover, 
the  details  of  the  reflected  waveform  must  take  into  account  reflections  from  all  of 
the  internal  stages,  in  order  to  match  in  complexity  the  experimental  conditions.  Thus, 
while  the  interpretation  of  the  data  shown  in  Figures  10  and  11  must  remain  qualitative, 
it  is  nevertheless  encouraging  to  find  that  the  overall  behavior  is  consistent  with  a 
model  which  predicts  that  an  incident  compression  reflects  from  the  inlet  as  an  expansion. 
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SECTION  5 

CONCLUDING  REMARKS 


As  pointed  out  above,  numerical  predictions  of  the  phenomena  measured 
in  experiments  with  real,  operating  engines  require  the  incorporation  of  many  effects, 
multiple  internal  reflections  being  one  of  the  most  important.  Unfortunately,  present- 
day  computer  codes  are  not  adequate  to  handle  this  problem.  Indeed,  the  present 
research  was  undertaken  in  an  effort  to  clarify  the  basic  physics  of  wave  transmission 
and  reflection  at  a  compressor  or  turbine  stage.  That  goal  has  been  achieved:  simple 
reflection  and  transmission  laws  have  been  derived,  and  are  qualitatively  in  accord  with 
the  measurements.  Moreover,  calculations  for  multi-stage  machines  can  be  done  by 
incorporating  the  actuator-disk  equations  as  a  means  of  linking  adjacent  stages. 

A  key  element  is  missing  from  current  understanding  of  the  physics, 
however,  namely  the  departure  of  the  mass-flow/pressure-ratio  relation  from  its  steady- 
state  value.  The  results  reported  here  have  all  been  based  on  the  assumption  that  the 
stagnation  pressure  ratio  can  be  specified  as  a  constant,  independent  of  the  mass  flow 
rate.  The  results  themselves  show  clearly  that  this  approximation  needs  to  be  improved, 
since  they  predict  mass  flow  changes  by  factors  on  the  order  of  1.5  during  the  passage 
of  a  wave.  Excursions  of  this  magnitude  exceed  the  limits  of  many  steady-state 
compressor  or  turbine  maps,  and  call  for  improved  modeling. 

Recent  studies  of  rotating  stall  and  surge  (12-13)  have  enjoyed  considerable 
success,  due  in  part  to  their  use  of  unsteady  stage-performance  characteristics,  in  which 
the  pressure  rise  Ap  is  taken  to  depend  on  both  the  flow  rate  ^  and  its  rate  of 
change  <p  : 

AP 

- r  ^  (<!>)-  ww 

f  ev 

The  time  constant  T  appearing  in  this  expression  is  related  to  the  time  required  to  fill 
the  plenum  volume  downstream  of  the  compressor  (the  burner  inlet,  for  example)  and 
so  this  relation  is  probably  not  applicable  to  the  present  problem,  where  the  mass  flow 
changes  occur  on  a  much  shorter  time  scale. 
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What  is  needed  is  a  combined  experimental  and  theoretical  study,  aimed 
at  identifying  the  instantaneous  mass  flow/pressure  rise  relation  for  a  single  stage. 
The  experimental  apparatus  to  accomplish  this  is  already  available  at  Calspan;  with 
only  slight  modifications,  it  would  be  possible  to  send  controlled  pressure  waves  through 
a  single  stage,  which  would  be  specially  instrumented  to  measure  the  (time-dependent) 
mass  flow  and  the  pressure  variations  in  the  transmitted  and  reflected  waves. 

In  parallel  with  these  experiments,  an  analytical  study  would  be  undertaken 
to  solve  for  the  detailed  modifications  of  the  flow  that  occur  when  a  plane  wave  passes 
through  a  stage.  The  work  required  amounts  to  a  calculation  of  shock  propagation 
through  a  pair  of  curved  channels,  whose  total  turning  angle  ii.  the  absolute  frame  of 
reference  is  zero.  The  analogous  problem  of  shock  propagation  through  em  area  change 
has  been  examined  by  Whitham,  Chisnell  and  Miles  (16-18)  and  has  led  to  extremely 
simple  relations  for  describing  flows  through  ducts  whose  cross-sectional  area  undergoes 
large  changes. 

The  clarification  of  the  pressure-ratio/mass  flow  relation  at  the  level  of 
a  single  stage  and  on  a  very  short  time  scale  is  the  key  ingredient  that  prevents  current 
state-of-the-art  codes  from  predicting  wave  passage  through  multistage  components. 
Once  it  is  known,  these  codes  can  be  revised  in  such  a  way  as  to  embrace  the  correct 
physics  at  the  basic  level. 

Knowledge  of  this  microscopic  relation  also  holds  the  key  to  the  relation 
between  these  short-duration  disturbances  and  the  longer-term  problems  of  surge  and 
rotating  stall.  Different  time  scales  are  involved  in  blast-generated  transients  from 
those  of  surge  and  stall,  and  it  is  necessary  to  understand  the  detailed  mechanics  of 
the  former  phenomena  in  order  to  clarify  how  these  disturbances  lead  into  the  latter  ones. 


SECTION  6 
RECOMMENDATIONS 


In  order  to  make  further  progress  beyond  the  qualitative  understanding 
described  in  this  report,  it  is  necessary  to  investigate  the  detailed  relation  between 
pressure  ratio  and  mass  flowrate  during  the  passage  of  a  pressure  transient.  This  will 
require  a  study  of  the  mechanics  of  pressure  propagation  through  a  rotating  stage,  on 
the  short  time  scale  given  by  the  stage  length  divided  by  the  speed  of  sound. 

The  basis  for  this  recommendation  is  as  follows:  there  are  two  principial 
mechanisms  by  which  the  ingestion  of  a  pressure  wave  may  damage  an  operating  engine: 

1)  the  overpressure  may  produce  structural  damage 

2)  a  condition  of  surge  or  rotating  stall  may  be  induced 

In  the  first  case,  the  results  of  this  report  support  the  general  results  of 
previous  measurements  at  Calspan:  the  transmission  coefficients  are  less  than  one, 
which  suggests  that  if  the  structure  ahead  of  the  compressor  face  can  sustain  the 
overpressure,  then  the  remaining  components  in  the  flow  path  will  presumably  also 
sustain  the  transient. 

However,  there  is  always  the  possibility  that  certain  geometries  (stage 
lengths,  for  example)  or  incident  wave  forms  might  lead  to  focusing  of  multiply  reflected 
waves,  or  that  massflow  amplification  might  produce  an  instantaneous  choking  that 
would  send  a  hammber  shock  upstream.  To  resolve  this  question,  it  is  necessary  to  do 
a  more  complete  job  of  linking  multiple  stages,  and  this  in  turn  requires  an  accurate 
model  of  the  time-dependent  pressure/massflow  relation  for  each  stage. 

For  the  second  mechanism  described  above,  the  recommended  research  is 
the  same.  Current  studies  of  rotating  stall  and  surge  make  extensive  use  of  a  time- 
dependent  pressure/massflow  relation  that  is  based  on  experiments  and  analyses  of  the 
much  longer  transients  associated  with  filling  of  the  plenum  chamber  downstream  of  a 
compressor.  In  order  to  close  the  gap  between  these  models  and  the  results  of  the 
present  report,  it  is  necessary  to  derive  the  comparable  relation  for  propagation  of  a 
pressure  wave  through  a  single  stage. 
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APPENDIX  A 

FORTRAN  PROGRAM  FOR  DOWNSTREAM-PROPAGATING  WAVES 


C  PROGRAM  WAVES  -  CALCULATES  THE  TRANSMISSION  AND  REFLECTION  Of 
C  PRESSURE  WAVES  AT  Af'l  ACTUATOR  DISK. 

C  WAVES  CAN  IW’JDLE  INCIDEfTT  EXP^^SIOT'IS  AND  SHOCKS 
0PEN(UNIT=2.F1LE='WAVE0UT'  .STATUS=-'NEW^ ) 

1  PRINT*, 'INPUT  STAGE  PRESSURE  RATIO  (99.  TO  STOP) .EMI .GAMMA' 
READ*, PR, EMI ,GM 
IF(PR.EQ.99.)  GO  TO  30 
GMINV  =  1 ./GM 
GMl  =  GM  -  1. 

GMllNV  =  1 ./GMl 
T8GM1  =  2. /GMl 
GB»11  =  GM/GMl 
GPl  =  GM  ♦  1  . 

GMIBG  =  GMl/GM 
GM1B2  =  GMl/2. 

GMIBTG  =  GMlBG/2. 

TAU  =  PR*«GM1BG 
WRITE<2.100)  PR. EMI .GM. TAU 

100  FORMAT <1  OX. 'SOLLmOr>l  FOR  PRESSURE  WAVES  INCIDDfT  ON  AN  ACTUATOR' 
«  '  DISK' ./lOX.'OF  STAGT^^T I  ON-PRESSURE  RATIO  =  '  , F5 . 1 ./I  OX . 

»  'INLET  MACH  NUMBER  ='.F5.2,'  SPECIFIC-HEAT  RATIO  ='.F5.2, 

«  /lOX.'STAGT'tATIOf'l-TEMPERATURE  RATIO  IS',F7.3./) 

FIND  CONDITIOf'IS  IN  REGIOT'I  4: 

EMISQ  =  EM1*EM1 
TOM  =  2. /GMl /EMISQ 
OPTOM  =  1 .  +  TOM 
U4BU1  =  1 . 

KOLf'IT  =  1 

5  COrTTINUE 
OLD  =  U4BU1 

ARG  =  (TAU*OPTOM  -  U4BU1  «U4BU1  )/'TOM 
IFCARG.GT.O.)  GO  TO  22 

PRINT*. 'NEGATIVE  ARGUMENT  IN  U4/UI  ITERATIONS' 

GO  TO  1 

22  U4BU1  =  ARG**(-GM1 INV) 

ERR  =  ABS(U4BU1-0LD) 

IFCERR.LT.l .E-07)  GO  TO  6 
KOUr>IT  =  KOUNT  +  1 

iF<Kour'nr.LE,30)  go  to  s 

PRirfr*.'U4/Ul  ITERATIOT'IS  FAILED  TO  COT'IVERGE 
GO  TO  1 

6  R4BR1  =  1 ./U4BU1 
P4BP1  *  R4BR1**GM 
T4BT1  *  P4BP1/R4BR1 
A4BA1  »  SQRT<T4BT1) 

TlBTOl  =  l./(l.  ♦  GM1B2»EM1SQ) 

PlBPOl  =  T1BT01**GBGM1 
EM4  =  U4BU1*EM1/A4BA1  35 

U4BA1  a  EM4*A4GA1 


oooo  oooo 


P4BP04  =  (1.  +  GM1B2*EM4«EM4)*«-GBGM1 

PRINT*. 'U4/A1  R4/R1  P4/P1  T4/T1  A4/A1  M4  P04/P01' 

PRINT*, U4BA1 ,R4BR1 ,P4BP1 .T4BT1 ,A4BA1 ,EM4,PR 
IREG  =  4 

WRITE<2,101)  IREG,U4BA1 ,R4BR1 ,P4BP1 ,T4BT1 ,A4BA1 ,EM4,PR 

101  F0RMAT< IX. 'REGION  U4/A1  R4/R1  P4/P1  T4/T1  A4/A1  M4' , 

*  '  P04/P01' ./I3.3X,7F7.4) 

CHOOSE  P2/P1.  AND  FIND  THE  CORRESPONDING  CONDITIONS  IN  REGIOT'I  2: 

PRU'IT*.' INPUT  P2/P1  <99.  TO  STOP)' 

READ*,P2BP1 

IF(P2BP1 .EQ.99.)  GO  TO  30 
IF(P2BP1  .GE.l .)  GO  TO  565 
INCIDENT  UAUE  IS  AN  EXPANS lOf'l: 

USIBAI  =  1 .  ♦  EMI 
R2BR1  =  P2BP1«*GMINV 
T2BT1  =  P2BP1/R2BR1 
A2BA1  =  SQRT<T2BT1) 

P02B01  =  1. 

U2BA1  =  TBGM1»<A2BA1-I .)  +  EMI 
EM2  =  U28A1/A2BA1 
US2BA1  =  A2BA1  *  U2BA1 
GO  TO  566 

565  USIBAI  =  EMI  ♦  SQRT< <GP1 *P2BP1 ♦GMl )/2 ./GM) 

R2BR1  =  <GP1*P2BP1  +  GMl )/(GMl »P2BP1  +  GPl) 

T2BT1  =  P2BP1/R2BR1 
A2BA1  »  SQRT<T2BT1) 

P02B01  =  <P2BP1**<-GM1IW))«(R2BR1«*GBGM1) 

EM2  =  US1BA1/A2BA1  -  SQRT< <GM1 «P2BP1 ♦GPl >/<2.«GM*P2BPl ) > 

U2BA1  =  EM2*A2BA1 

566  PRir’rr*.'USl/Al  U2/A1  R2/R1  P2/P1  T2.'T1  A2/A1  M2  P02/P01' 

PRINT*. USIBAI .U2BA1 .R2BR1 ,P2BP1 ,T2BT1 ,A2BA1 ,EM2,P02B01 

IREG  =  2 

WRITE':2.102)  IREG. USIBAI  ,U2BA1  .R2BR1  .P2BP1  .T2BT1  ,A2BA1  ,EM2, 

*  P02B01 

102  F0RMAT</1X. 'REGION  USl/Al  U2/A1  R2/R1  P2/P1  T2/T1  A2/Ar , 

*  '  M2  P02/P01' ,/I3,3X.8F7.4) 

IF<US2BA1 .NE.O.)  WRITE<2.115)  US2BA1 

115  FORMAT <3X,'US2/A1  =',F7.4) 

NOW  ITERATE  Or>l  P5/P4  Al-JD  P3/P2  UNTIL  PRESSURES  AND  PARTICLE 
VELOCITIES  IN  REGIONS  5  AND  6  ARE  MATCHED 

ITER  =  1 
KASE  a  1 
P5BP4  =  1 . 

P54A  a  1  . 

P54B  a  P2BP1 
P3BP2  a  1 . 

P32A  a  1  . 

P32B  a  P2BP1 
10  COf'ITINUE 
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IF(P2BP1 .6T.1 .)  60  TO  301 
C  FOR  THE  CASE  OF  AN  INCIDENT  EXPANSION; 

R5BR4  =  P5BP4**GMINW 
T5BT4  *  P5BP4/R5BR4 
A5BA4  =  SQRT(T5BT4) 

US5BA4  «  U4BA1  *  A4BA1 
P05B04  =  1  . 

U5BA4  =  TBGM1*(A5BA4-1  .)  ♦  EM4 
EM5  =  U5BA4/A5BA4 
US56A5  -  US5BA4/A5BA4 
GO  TO  302 

C  FOR  THE  CASE  OF  AN  INCIDENT  SHOCK; 

301  US5BA4  =  EM4  ♦  SQRT< <GP1«P5BP4  ♦  GMl)/2./GM) 

R5BR4  =  <6P1«P5BP4+GM1)/<GM1*P5BP4*GP1) 

T5BT4  =  P5BP4/R5BR4 
A5BA4  =  SQRTa5BT4) 

P05B04  =  <P5BP4**<-GM1INU))»<R5BR4*»GBGM1) 

US5BA5  =  US5BA4/A56A4 

EM5  =  US5BA5  -  SQRT< <GM1»P5BP4  ♦  GPl )/2./GM/P5BP4) 

302  USBAl  s  EM5»A5BA4»A46A1 
T5BT1  =  T5BT4«T4BT1 
A5BA1  =  SQRT<T5BT1) 

USSBAl  -  USSBA5«A5BAi 
R5BR1  =  R5BR4*R4BR1 
P5BPI  =  P5BP4«P4BP1 
P05B01  =  P05B04*PR 

PRINT#/ US5/A1  U5/A1  R5/R1  P5/P1  T5/T1  A5/A1  EM5  P05/P01' 

PRINT#, USSBAl , USBAl .RSBRl .PSBPl ,T5BT1 .ASBAl ,EM5,P0SB01 
IF(KASE.EQ.2)  60  TO  11 
XSA  *  PSBPl 
YSA  a  USBAl 
GO  TO  63 
11  X5B  =  PSBPl 
Y5B  =  USBAl 

CHOOSE  P3/P2,  AND  THE  CORRESPOPl'JDING  C01>IDITI0NS  IN  REGIONS  3: 

63  IF<P3BP2.GT.0.)  GO  TO  563 

PRIf>(T#/P3/P2  NEGATIk^E  -  PRESS  RETURT'I  TO  COTfriNUE' 

PAUSE 

WRITE':2.110) 

no  FORMAT (/5X/P3/P2  NEGATIVE  -  GOING  TO  NEXT  CASE',//) 

GO  TO  1 

563  IF(P3BP2.LT.l .)  60  TO  60 

US3BA2  =  -  EM2  ♦  SORT (< GPl »P3BP2+GM1 )/2 ./GM) 

R3BR2  a  <GP1«P3BP2+GM1)/<6M1«P3BP2*GP1) 

T3BT2  a  P3BP2/R3BR2 
A3BA2  a  SQRT<T3BT2) 

P03B02  a  <P3BP2*»-GM1INV)»<R3BR2»#GBGM1> 

US36A3  a  US3BA2/A36A2 

EM3  a  -  US3BA3  ♦  SQRT<<GM1«P3BP2+GP1 )/2./GM/P3BP2) 

T3BT1  a  T3BT2#T2BT1 
A3BA1  a  SQRT<T3BT1) 

US3BA1  a  US3BA3#A3GA1 
R3BR1  a  R3BR2»R2BR1 
P3BP1  a  P38P2#P2BP1 
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P03B01  =  P03B02«P02B01 
U3BA1  »  EM3«A3BA1 
60  TO  66 

60  A3BA2  =  P3BP2«*GM1BTG 
A3BA1  =  A36A2»A2BA1 
T3BT1  *  A3BA1*A3BA1 
U3BA1  =  TBGH1«(A2BA1-A3BA1)  *  U2BA1 
P3BP1  =  P3BP2»P2BP1 
US3BA1  s  -A2BA1 
R3BR2  =  P3BP2**GMlNy 
R3BR1  =  R3BR2»R2BR1 

EM3  =  U3BA1/A3BA1  ? 

P03B01  =  P02B01 
66  cor^rriNUE 

PRir>rr».'US3/Al  U3/A!  R3/R1  P3/P1  T3/TS  A3BA1  EM3  P03/PC 

PRINT*, US3BA1 ,U3BA1 .R3BR1 .P3BP1 ,T3BT1 ,A3BA1 ,EM3,P03B01 

NOU  DO  THE  TRP#'JSITIOr<l  FROM  REGION  3  TO  REGION  6: 

P06B01  =  PR*P03B01 
EM3SQ  =  EM3*EM3 
TOM  =  2./GM1/EM3SQ 
OPTOM  =  1 .  ♦  TOM 
U<6BU3  =  1  . 

KOUNT  =  1 

55  COfTTINUE 
OLD  «  U<&BU3 

ARG  =  <(TAU*OPTOM  -  U6BU3*U6BU3)/T0M) 

IF<ARG.6T.O.)  GO  TO  52 

PRINT*. 'NEGATIVE  AR6UMEr>fr  IN  U6/U3  ITERATIONS' 

GO  TO  1 

52  U6BU3  =  AR6»*';-6MlINv') 

ERR  =  ABS(U6BU3-0LD) 

IF(ERR.LT.l .E-07)  GO  TO  56 
KOWT  =  KOUffT  ♦  I 
IF<KOUr'IT.LE.30)  GO  TO  55 

PRir>rr*,'U6/u3  iterations  failed  to  converge' 

GO  TO  1 

56  R6BR3  =  1 ./U6BU3 
P6BP3  =  R6BR3**GM 

T6BT3  =  P6BP3/R6BR3 
A6BA3  =  SQRT(T6BT3) 

EM6  =  U6BU3*EM3/A6BA3 
U66A1  »  U6BU3*U36A1 
T6BT1  =  T6BT3*T3BT1 
A6BA1  =  SQRT(T6BT1) 

R6BR1  =  R6BR3*R3BR1 
P6BP1  =  P6BP3*P3BP1 

PRINT*, 'U6/A1  R6/R1  P6/P1  T6/T1  A6/AI  M6  P06/P0r 

PRINT*, U6BA1 ,R6BR1 .P6BP1 ,T6BT1 ,A6BA1 ,EM6,P06B01 
IF<KASE.EQ.2)  60  TO  61 
X6A  *  P6BP1 
Y6A  =  U6BA1 
DX  =  X5A  -  X6A 
DY  =  Y5A  -  Y6A 
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DIFA  =  DX»DX  +  DY*DY 
KASE  =  2 
P5BP4  »  P54B 
P3BP2  =  P32B 
GO  TO  10 
61  X6B  =  P6BP1 
Y6B  s  U6BA1 

NOW  FIND  THE  NEXT  GUESS,  AND  ITERATE  ONE  MORE  TIME: 

DX  =  X5B  -  X6B 
DY  =  Y5B  -  Y6B 
DIFB  =  DX*DX  ♦  DY«OY 
DIF  =  SQRT(DIFB> 

IFCDIF.LT.l .E-04)  GO  TO  70 
ITER  =  ITER  ♦  1 
IF<ITER.LE.30)  GO  TO  564 

PRINT»,'ITERATia'JS  Of>l  P3/P2  AND  P5/P4  NOT  COhWERGING  -  ' 
PRINT*, 'PRESS  RETURN  TO  RESTART' 

PAUSE 

URITE<2.111) 

111  F0RMAT(2X,'ITERATI0r>IS  Of-I  P3./P2  AND  P5/P4  NOT  CONVERGENT',/) 

GO  TO  1 

564  SLP5  =  <Y5B-Y5A)/(X5B-X5A) 

B5  =  Y5A  -  SLP5*X5A 
SLP6  =  (Y6B-Y6A)/<X6B-X6A) 

B6  *  Y6A  -  SLP6»X6A 

Xf-IEW  =  -<B5-B6)/<SLP5-SLP6) 

P54NEU  =  P54A  ♦  <P54B-P54A)*(XNEW-X5A)/<X5B-X5A) 

P32NEW  =  P32A  +  <P32B-P32A)*<;XNEW-X6A)/(X6B-X6A) 

PRINT*, 'ITER  P5/P4  P3/P2  DIF' 

PRINT*,  ITER,  P54NEW.P32f'IEW,  DIF 
IF<U6BA: .LT.O.)  PRINT*, 'NOTE  -  U6  IS  NEGATIVE' 

PAUSE 

IF<DIFA.LT.DIFB)  GO  TO  71 
X5A  =  X5B 
X6A  =  X6B 
Y5A  =  Y5B 
Y6A  =  Y6B 
P54A  =  P54B 
P32A  =  P32B 
DIFA  =  DIF*DIF 
71  P54B  =  P54NEW 
P32B  =  P32f^EW 
P5BP4  »  P54B 
P3BP2  =  P32B 
GO  TO  10 
70  COI4TINUE 

WRITE(2,112)  P5BP4 

112  F0RMAT(/10X,'P5/P4  =',F7.4) 

IREG  «  5 

WRITE(2,104)  1REG,US5BA1 ,U5BA1 ,R5BR1 ,P5BP1  ,T5BT1  ,A5BA1 ,EM5, 
^  PO  5B0 1  iQ 


104  FORMAK IX /REGION  US5/A1  U5/A1  R5/R1  P5/P1  T5/T1  AS/Al", 

«  '  EM5  P05/P01' ,./I3,3X,8F7.4) 

IREG  =  3 

URITE-:2 .105)  P3BP2,IREG,US3BA1  ,U3BA1  ,R3BR1  ,P3BP1  ,T3BT1  ,A3BA1 , 

*  EM3,P03B01 

105  FORMAT(/10X,'P3/P2  =',F7.4, 

*  /IX, 'REGION  US3/A1  U3/A1  R3/R1  P3/P1  T3/T1  A3/A1  EM3' 

»  '  POG/POl' ,/I3,3X,8F7.4) 

IREG  =  6 

UIRITE(2,106)  IREG.U<&BA1  ,R^BR1  ,P6BP1  ,T6BT1  ,A6BA1  ,EN6.P06601 
104  F0RMAT</1X/REGI0N  U4/A1  R4/R1  P4/P1  T6/T1  A6/A1  M4' , 

*  '  P06/P01' ./I3,3X,7F7.4) 

EMDOTl  =  EMI 

EMD0T3  =  U3BAUR3BR1 
URITE(2,108)  EMDOTl  .E.MD0T3 
108  FORMAT (1  OX, 'MASSFLOU  IN  REGIONS  1  AND  4  =',F7.4, 

»  /lOX.'MASS  FLOW  IN  REGIOr<IS  3,5,  AND  6  »',F7.4,//) 

GO  TO  1 
30  COf^TINUE 

CL0SE<Ur'IIT=2) 

STOP 
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APPENDIX  B 

FORTRAN  PROGRAM  FOR  UPSTREAM-PROPAGATING  WAVES 

C  PROGRAM  WAVEC  -  CALCULATES  THE  TRANSMISSION  REFLECTIOf^  OF 
C  PRESSURE  UAL>ES  AT  AN  ACTUATOR  DISK. 

C  UIAVEC  H^DLES  EXP^SIONS  AND  SHOCKS  INCIDENT  FROM  THE  DOWNSTREAM 
C  SIDE  OF  THE  ACTUATOR  DISK 

OPEN<UNIT=2,FlLE»'WAL»ECOUT"  ,STATUS='NEW' ) 

1  PRINT*, 'INPUT  STAGE  PRESSURE  RATIO  (99.  TO  STOP ) , EM I , GFW1A ' 
READ*, PR. EMI ,GM 
IF(PR.EQ.99.)  GO  TO  30 
GMINV  =  1 ./GM 
GMl  =  GM  -  1. 

GMIIW  =  1  ./GMl 
TBGM!  »  2. /GMl 
GBGMl  =  GM/'GMl 
GPl  =  GM  +  1. 

GM1B6  =  GMl/GM 
GM1B2  =  GMl/2. 

GMIBTG  =  GMlBG/2. 

TAU  *  PR*»GM1BG 
WRITE(2,100)  PR.EMl .GM.TAU 

100  FORMAT (1  OX. 'SOLUTION  FOR  PRESSURE  WAVES  INCIDENT  FROM  THE', 

«  '  DOWNSTREAM  SIDE  OF  AN  ACTUATOR', 

*  '  DISK' ,/10X,'OF  STAGNATION-PRESSURE  RATIO  =' , F5 . 1  ,/l OX , 

«  'INLET  MACH  NUMBER  =',F5.2,'  SPECIFIC-HEAT  RATIO  *',F5.2. 

*  / 1  OX, 'STAGT^IAT  ION-TEMPERATURE  RATIO  IS',F7.3,/) 

FIND  COf'JOITIONS  IN  REGIOf'l  4; 

EMISQ  =  EM1»EM1 
TOM  =  2./GM1/EM1SQ 
OPTOM  *  1 .  +  TOM 
U48U1  =  1 . 

KOUNT  =  1 

5  COf'iTINUE 
OLD  =  U4BU1 

AR6  =  <TAU*OPTOM  -  U4BU1 «U4BU1 )/TOM 
IF<ARG.GT.O.)  GO  TO  22 

PRINT*, 'NEGATIVE  ARGUMENT  IN  U4/U1  ITERATIONS' 

GO  TO  1 

22  U4BU1  »  ARG**<-GM1INV) 

ERR  =  ABS<U4BU1-0LD) 

IF<ERR.LT.l .E-07)  GO  TO  6 
KOUNT  *  Kourrr  ♦  i 
IF(KOLfTT.LE.30)  GO  TO  5 

PRINT*. 'U4/U1  ITERATIONS  FAILED  TO  COT^VERGE ' 

GO  TO  1 

6  R4BR1  =  1./U4BU1 
P4BP1  »  R4BR1**GM 
T4BT1  =  P4BP1/R4BR1 
A4BA1  =  SQRT(T4BT1) 
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TlBTOl  *  +  QM1B2*EM1SQ) 

PlBPOl  *  TlBTOHn»GBGMl 
EM4  s  U4BUUEMI/A4BA1 
U4BA1  »  EM4*A4BA1 

P4BP04  =(1.4  GM1B2*EM4»EM4)*#-GBGM1 
PRINT*, 'U4/A1  R4/R1  P4/P1  T4/T1  A4/A1  M4  P04/P0r 

PRINT*, U4BA1 ,R4BR1 ,P4BP1 ,T4BT1 ,A4BA1 ,EM4,PR 
I  REG  »  4 

WRITE(2,101)  IREG,U4BA1 ,R4BR1 .P4BP1 ,T4BT1 ,A4BA1 ,EM4,PR 

101  FORPWT( IX, 'REGION  U4/A1  R4/R1  P4/P1  T4/T1  A4/A1  M4', 

«  '  P04/P0r  ,/I3,3X,7F7.4) 

CHOOSE  P5/P4,  AND  FIND  THE  CORRESPONDING  CONDITIONS  IN  REGIW'I  5; 

PRINT*, 'INPUT  P5/P4  <99.  TO  STOP>' 

READ*,P5BP4 

IF<P5BP4.EQ.99.)  GO  TO  30 
IF<P5BP4.GE.l .)  GO  TO  565 
C  INCIDENT  WAOE  IS  Af-l  EXPANSIOf^l: 

US5BA4  =  014  -  1 . 

P5BP1  =  P5BP4*P4BP1 
R5BR4  =  P5BP4*»GMIW 
R5BR1  =  R5BR4*R4BR1 
T5BT4  =  P5BP4/R5BR4 
T5BT1  =  T5BT4*T4BT1 
A5BA4  =  SQRT(T5BT4) 

A5BA1  s  A5BA4*A4BA1 
P05B04  =  1  . 

P05B01  «  PR 

U5BA4  =  TBGM1*<1 .-A5BA4)  ♦  EM4 
U5BA1  ==  U5BA4*A4BA1 
EM5  =  U5BA4/A5BA4 
GO  TO  566 

C  INCIDENT  WAVE  IS  A  SHOCK: 

565  US5BA4  =  EM4  -  SORT-:  <GP1*P5BP4+GMI  )/2./GM) 

R5BR4  =  <GP1*P5BP4  ♦  GMl )/<GMl *P5BP4  +  GPl) 

R5BR1  =  R5BR4*R4BR1 
P5BP1  =  P5BP4*P4BP1 
T5BT4  =  P5BP4/R5BR4 
T5BT1  =  T5BT4*T4BT1 
A5BA4  =  SQRT<T5BT4) 

A5BA1  =  A5BA4»A4BA1 

P05B04  =  (P5BP4*»(-GM1 INV))*(R5BR4*»GBGM1 ) 

P05B01  =  P05B04*PR 
US5BA5  =  US5BA4/A5BA4 

EM5  =  US5BA5  +  SQRT< <GM1*P5BP4+GP1 )/< 2  .*GM*P5BP4) ) 

U5BA1  =  EM5»A5BA1 

566  PRINT*. 'US5/A4  U5/A1  R5/R1  P5/P1  T5/T1  A5/A1  M5  P05/P01' 

PRINT*,US5BA4,U5BA1 .R5BR1 ,P5BP1 ,T5BT1 ,A5BA1 ,EM5,P05B01 

IREG  =  5 

WRITE<2,102)  IREG.US5BA4,U5BA1 ,R5BR1 ,P5BP1 ,T5BT1 .A5BA1 ,EM5, 

*  P05B01 

102  F0RMAT(/1X. 'REGION  US5/A4  U5/A1  R5/R1  P5/P1  T5/T1  A5/A1', 

*  '  M5  P05/P01' ,/I3,3X,8F7.4) 
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NOW  ITERATE  ON  P2/P1  P6/P5  UNTIL  PRESSURES  AND  PARTICLE 

VELOCITIES  IN  REGIONS  6  7  ARE  hWTCHED 

ITER  s  I 
KASE  »  I 
P2BP1  »  1 . 

P21A  =  1  . 

P21B  »  P5BP4 
P4BP5  *  1 . 

P65A  =  1  . 

P65B  =  P5BP4 
JO  CONTINUE 

IF<P5BP4.GT.l .)  GO  TO  301 
C  FOR  THE  CASE  OF  Pt^  INCIDENT  EXPFWSION: 

R2BR1  =  P2BP1»»6M1NV 
T2BT1  =  P2BP1/R2BR1 
A2BA1  =  S0RT<T2BT1) 

US2BA1  =  EMI  -  1. 

P02B01  =  1 . 

U2BA1  =  TBGMUKl  .-A2BA1)  ♦  EMI 
EM2  =  U2BA1/A2BA1 
LIS2BA2  =  US2BA1/A2BA1 
60  TO  302 

C  FOR  THE  CASE  OF  AN  INCIDENT  SHOCK; 

301  US2BA1  =  EMI  -  SQRT< < GP1*P2BP1  +  QMl)/2./GM) 

R2BR1  =  <GP1»P2BP1*6M1)/<GM1*P2BP1+GP1> 

T2BT1  =  P2BP1/R2BR1 
A2BA1  =  SQRT(T2BT1) 

P02B01  =  (P2BP1«*(-GM1INV))*(R2BR1»»6BGM1) 

US2BA2  »  US2BA1/A2BA1 

EM2  =  US2BA2  ♦  SQRT< <6M1«P2BP1  ♦  6PI)/2 ./GM'P2BP1 ) 

302  U2BA1  =  EM2»A2BA1 

PRINT*, 'US2/A1  U2/A1  R2/R1  P2/P1  T2/TI  A2/A1  EM2 

PRINT*, US2BA1 ,U2BAI ,R2BRI ,P2BP1 .T2BT1 ,A2BA1 ,EM2.P02B01 

C  NOW  DO  THE  TRANSITIOTV  FROM  REGION  2  TO  REGION  7; 

C 

P07B01  =  PR»P02B01 
EM2SQ  =  EM2*EM2 
TOM  a  2./GM1/EM2SQ 
OPTOM  a  1 .  ♦  TOM 
U7BU2  a  I  . 

KOll'IT  a  1 
55  COf'TriNUE 
OLD  =  U7BU2 

ARG  a  (<TAU«0PT0M  -  U7BU2«U7BU2)/T0M) 

IF(ARG.GT.O.)  GO  TO  52 

PRINT*, 'NEGATIVE  ARGUMEITT  IN  U7/U2  ITERATIONS' 

GO  TO  1 

52  U7BU2  a  AR6««(-GM1INL 
ERR  a  ABS<U7BU2-0LD) 

IF<ERR.LT.l .E-07)  GO  TO  56 
KOUT'IT  a  kOUNT  +  1 
IF(KOUNT.LE.30)  GO  TO  55 

PRINT*, 'U7/U2  ITERATIOr^S  FAILED  TO  COhWERGE' 

GO  TO  1 


P02/P0] ' 
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56  R7BR2  *  1  ./U7BU2 
P7BP2  »  R7BR2*«GM 
T7BT2  =  P7BP2/R7BR2 
A7BA2  =  SQRT<T7BT2) 

EM7  =  U7BU2«EM2/A7BA2 
U7BA1  =  U7BU2*U2BA1 
T7BT1  =  T7BT2«T2BT1 
A7BA1  =  SQRT(T7BT1) 

R7BR1  =  R7BR2»R2BR1 
P7BP1  »  P7BP2*P2BP1 

PRINT*. 'U7/A1  R7/R1  P7/P1  17/11  A7/A1  M7  P07/P0r 

PRINT*, U7BA1 ,R7BR1 ,P7BP1 ,T7BT1 ,A7BA1 .EM7,P07B01 
IF<i<ASE.EQ.2)  60  TO  11 
X7A  s  P7BP1 
Y7A  =  U7BA1 
X2A  =  P2BP1 
Y2A  =  U2BA1 
GO  TO  63 

11  X7B  =  P7BP1 
Y7B  =  U7BA1 
X2B  =  P2BP1 
Y2B  *  U2BA1 

CHOOSE  P6/P5 

63  IF<P6BP5.GT,Q.)  60  TO  563 

PRIf4T*,'P6/P5  NEGATIVE  -  PRESS  RETURN  TO  CONTINUE' 
PAUSE 

URITE<2.110) 

no  FORMAT </5X.'P6/P5  NEGATIVE  -  GOING  TO  NEXT  CASE',//) 

GO  TO  1 

563  IF<P6BP5.LT.l .)  GO  TO  60 

US68A5  =  EM5  +  SQRT< (GPl *P6BP5*GM1 )/2./GM) 

R6BR5  =  <GP1«P6BP5+6M1)/(GM1*P6BP5+GP1) 

T6BT5  =  P6BP5/R6BR5 
A6BA5  =  SQRT<T6BT5) 

P06B05  =  ':P6BP5**-GN1IW)*<R6BR5*»GBGM1) 

US6BA6  »  US6BA5/A66A5 

EM6  =  US6BA6  -  SQRT( (6M1»P6BP5+GP1 )/2 ./GM/P6BP5) 
T6BT1  =  T6BT5«T5BT1 
A6BA1  =  SQRT<T6BT1) 

US6BA1  -  US66A6«A6BA1 
R6BR1  =  R6BR5»R5BR1 
P6BP1  «  P6BP5*P5BP1 
P06B01  *  P06B05*P05B01 
U6BA1  s  EM6*A66A1 
60  TO  66 

60  A6BA5  =  P6BP5»*GM1BTG 
A6BA1  >  A6BA5»AS6A1 
T6BT1  »  A6BA1*A6BA1 
U6BA1  *  TB6M1»<A6BA1-A5BA1)  ♦  U5BA1 
P6BP1  *  P6BP5*P5BP1 
US6BA1  s  -A56A1 
R6BR5  =  P6BP5**GMINV 
R6BR1  >  R6BR5*R5BR1 
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EM6  s  U6BA1/A6BA1 
P06B01  =  P05B01 
66  CONTINUE 

PRINT»,'US<6/A1  U6/A1  R<5/R1  P6/Pi  16/11  A6BA1  EM6  P06/P0r 

PRINT*, US6BA1  ,U6BA1 ,R6BR1 ,P6BP1 ,T<5BT1 ,A6BA1 ,EM6,P06B01 

IF(KASE.EQ.2)  GO  TO  61 
X6A  =  P6BP1 
Y6A  s  U6BA1 
DX  =  X6A  -  X7A 
DY  =  Y6A  -  Y7A 
DIFA  =  DX»DX  +  DY*DY 
KASE  =  2 
P2BP1  =  P21B 
P6BP5  =  P65B 
GO  TO  10 
61  X6B  =  P6BP1 
Y6B  =  U6BA1 

NOW  FIND  THE  NEXT  GUESS.  AND  ITERATE  OI'IE  MORE  TIME: 

DX  =  X6B  -  X7B 
DY  =  Y6B  -  Y7B 
DIFB  =  DX»DX  +  DY*DY 
DIF  =  SQRT(DIFB) 

IF<D1F.LT.1  .E-0<;)  GO  TO  70 
ITER  =  ITER  ♦  1 
IF(ITER.LE.30)  G'  TO  564 

PRINT*,  ■'ITERATIOT'IS  ON  P6/P5  Af^lD  P2/P1  NOT  CONVERGING  -  ' 

PRINT*, 'PRESS  RETURf'l  TO  RESTART' 

PAUSE 

WRITE';  2,  111) 

111  F0RMAT(2X.' ITERATIONS  Of'l  P6/P5  AND  P2/P1  NOT  COr^lVERGENT' ./) 

GO  TO  1 

564  SLP6  =  <Y6B-Y6A)/<X6B-X6A) 

B6  =  Y6A  -  SLP6«X6A 
SLP7  =  (Y7B-Y7A)/(X7B-X7A) 

B7  =  Y7A  -  SLP7*X7A 

XNEU  =  -':B6-B7)/<SLP6-SLP7) 

P71NEW  =  Xt4EU 

P21NEW  =  P21A  +  <XNEW-X7A)*(P21B-P21A)/<X7B-X7A) 

P65NEW  =  XNEU/P5BP1 
PRINT*, 'ITER  P2/P1  P6/P5  DIF' 

PRINT*, ITER, P21NEU,P65NEW, DIF 
IF(U6BA1  .LT.O.)  PRINT*, 'NOTE  -  U6  IS  NEGATIVE' 

PAUSE 

IFCDIFA.LT.DIFB)  GO  TO  71 
X6A  =  X6B 
X7A  =  X7B 
Y6A  =  Y6B 
Y7A  =  Y7B 
P21A  =  P21B 
P65A  =  P65B 
DIFA  =  DIF*DIF 
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71  P21B  =  P21NEU 
P<45B  =  P<65NEUI 
P2BP1  =  P21B 
P6BP5  =  P45B 
GO  TO  10 
70  CONTINUE 

WRITE<2.112)  P2BP1 
112  FORMAT(/10X/P2/P1  =',F7.4) 

I  REG  =  2 

U)R]TE<2.104)  IRE6,US2BA1  ,U2BA1  .R2BR1  ,P2BP1  ,T2BT1  ,A2BA1  ,EM2, 

*  P02B01 

104  FORMAT'MX.'REGIOl'l  US2/A1  U2/A1  R2/R1  P2/P1  T2./T1  A2/Ar, 

*  '  EM2  P02/P01' ,/I3,3X,8F7.4) 

IREG  =  6 

WRITE<2,105)  P6BP5,lREG.US<iBAl  ,U6BA1  ,R(SBR1  ,P<4BP1  ,T6BT1  ,A6BA1  , 

*  EM6,P0<iB01 

105  FORMAT  </l  OX. 'P(i/P5  =',F7.4, 

*  /lX.'REGIOt'l  US6/A1  U6/A1  R6/R1  P6/P1  T6/T1  A<S/A1  EM6' 

*  '  P06/P0r  ,/I3,3X,8F7.4) 

IREG  =  6 

UJRITE<2,106)  1RE6.U6BA1  ,R6BR1  ,P6BP1  ,T(SBT1  ,A6BA1  ,EM6,P06B01 

106  FORMAT (/IX.'REGIOf'l  U6/A1  R6/R1  P6/P1  T6/T1  A6/A1  M6', 

«  •'  P06/P0r  ,/I3,3X,7F7.4) 

EMDOTl  =  EMI 
EMD0T2  =  U2BA1»R2BR1 
WRITE': 2, 108)  EMDOTl  ,EMD0T2 
108  FORMAT  (1  OX, 'MASSFLOW  IN  RE6IOr>IS  1  AND  4  s',F7.4, 

*  /lOX/MASS  flow  in  REGIONS  2,7,  AND  6  =',F7.4,//) 

GO  TO  1 

30  CW'ITINUE 

CLOSED UNIT=2) 

STOP 

END 
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